##[PUR]=group
##wd=folder
##lu=raster
##zo=raster
##loolu=file
##loozo=file
##passfilenames
##showplots

library(raster)
library(rasterVis)
library(reshape2)
library(plyr)
library(lattice)
library(latticeExtra)
library(RColorBrewer)
library(hexbin)
library(grid)
library(ggplot2)
library(R2wd)
library(foreign)

#set project properties
setwd(wd)

#load datasets (land use t1, land use t2, zone)
landuse1 <- raster(lu)
zone <- raster(zo)
Lookup_LU<-read.table(loolu, header=TRUE,sep=",")
Lookup_Zone<-read.table(loozo, header=TRUE,sep=",")

# set raster attribute table (RAT)
landuse1<-ratify(landuse1, filename='landuse1.grd',count=TRUE,overwrite=TRUE)
zone<-ratify(zone, filename='ratify.grd',count=TRUE,overwrite=TRUE)

#create land use change database
area_lc1<-as.data.frame(levels(landuse1))
area_zone<-as.data.frame(levels(zone))
area_x<-sum(area_zone$COUNT)
levels(landuse1)<-merge((levels(landuse1)),Lookup_LU,by="ID")
levels(zone) <- merge(area_zone,Lookup_Zone,by="ID")
area_lc1<-as.data.frame(levels(landuse1))
area_zone<-as.data.frame(levels(zone))
colnames(area_lc1)[2] = "COUNT_LC1"
colnames(area_lc1)[1] = "ID_CLASS"
colnames(area_zone)[2] = "COUNT_ZONE"
colnames(area_zone)[1] = "ID_ZONE"
cross <- as.data.frame(crosstab((stack(landuse1,zone))))
colnames(cross)[1] = "ID_CLASS"
colnames(cross)[2] = "ID_ZONE"
colnames(cross)[3] = "COUNT_CROSS"
cross<-merge(cross,area_lc1, by="ID_CLASS")
cross<-merge(cross,area_zone, by="ID_ZONE")
cross$lprop_lclc<-cross$COUNT_CROSS/cross$COUNT_LC1    #calculate proportion of land cover by total for each land cover class at landscape level
cross$lprop_lczo<-cross$COUNT_CROSS/cross$COUNT_ZONE  #calculate proportion of landcover by total area zone
cross$lprop_lc_l<-cross$COUNT_CROSS/area_x #calculate proportion for landscape area

#analyze data
cross.melt <- melt(data = cross, id.vars=c('Zone','CLASS'), measure.vars=c('COUNT_CROSS'))
prop<- dcast(data = cross.melt, formula = Zone ~ CLASS, fun.aggregate = sum)
plot2<-ggplot(data=cross.melt)+aes(CLASS,value,fill=Zone)+geom_bar(stat="identity",position="stack")+theme(axis.text.x= element_text(angle=90,hjust=1))
plot2

#write report
wdGet()
wdWrite("This document is automatically generated by LUMENS")
wdSection("Land Use Planning for Multiple Environmental Services")
wdTitle("PUR:Planning Unit Reconciliation",label="R2wd")
wdWrite("PUR: Planning Unit Reconciliation adalah sub-modul bagian dari LUMENS (Land Use planning for Multiple ENvironmental Services) yang berfungsi untuk melakukan rekonsiliasi data-data unit perencanaan dari sebuah daerah",paragraph = TRUE)
wdTable(area_zone, caption = "Planning Unit Area", caption.pos="above", autoformat=2)
wdTable(area_lc1, caption = "Land USe System Area", caption.pos="above", autoformat=2)
wdPlot(plot2,caption="Proportion of Land USe System in Planning Unit Classes",height = 6.08, width = 6.08,pointsize = 6,paragraph = TRUE)
wdTable(prop, caption = "Proportion of Land USe System", caption.pos="above", autoformat=2)
